/**
 * AdaptiveArray neighbour iterators.
 *
 * \file aa_neighbour.hpp
 * \author tomby
 * \date 04-01-2009
 **/

/*  Copyright (C) 2009  Tomas 'tomby' Bily

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>. */

#ifndef __DBC_AA_NEIGHBOUR__
#define __DBC_AA_NEIGHBOUR__

#define PMOD1(X,Y) (((X) & ((1 << (Y)) - 1)) ? 1 : 0)

namespace dbc {

  /**
   * Neighbour array iterator for blitz::Array
   * Iteration is done only on axes
   *
   * example:
   * + - center site
   * @ - neighbour sites of +
   *
   *        @
   *        @ 
   *      @@+@@
   *        @
   *        @
   */

  template<typename type, int dim>
  class ArrayIterator
  {
  public:
    typedef type Type;
    typedef blitz::Array<type, dim> ArrayType;
    typedef blitz::TinyVector<int, dim> IndexType;
    typedef blitz::RectDomain<dim> RectDomType;
    typedef RectDomIterator<dim> RectDomIter;
    enum { Dim = dim };

  protected:
    /// Array where we iterate
    ArrayType*	array;
    /// position of center site - iteration os around this site
    RectDomType	center;
    /// current neighbour 
    RectDomIter pos;
    /// false - laft side, true - right side
    bool	sign;
    /// active dimension
    int		active_idx;
    /// Boundary of iteration
    RectDomType	radius;
    /// true - if iteration was clipped by boundary of array, false - otherwise 
    bool clipped;

  public:
    ArrayIterator ();

    /**
     * Constructor
     *
     * \param iarray - array where we will iterate
     * \param icenter - center site 
     * \param radius - boundary of neighbourhood iteration (min,max)
     * \param end - indicate if iterator is set as finished
     */
    ArrayIterator (ArrayType* iarray, RectDomType& icenter, RectDomType& radius, bool end = false);

    /** 
     * Constructor for iteration in one dimension
     * 
     * \param iarray - array where we will iterate
     * \param icenter - center site 
     * \param radius - boundary of neighbourhood iteration (min,max)
     * \param a_idx - index of dimension where to iterate (0 .. dim - 1)
     */
    ArrayIterator (ArrayType* iarray, RectDomType& icenter, RectDomType& radius, int a_idx);

    /// Copy constructor
    ArrayIterator (const ArrayIterator& it);

    ArrayIterator& operator++ ();
    ArrayIterator operator++ (int);
    Type& operator* ();
    bool operator== (const ArrayIterator&);
    bool operator!= (const ArrayIterator&);

    /// operator ++ in one dimension
    ArrayIterator& pp ();
    IndexType position ();
    IndexType& rel_position ();
    IndexType& center_position ();
    Type& center_element ();
    int active_dim ();
    bool valid_p ();
    void invalidate ();
    RectDomType index_boundary ();
    void set_active_idx (int idx);
    ArrayType* data ();

    inline bool end_p ();
    inline bool clipped_p ();

    void print_dbg ();
  protected:
    inline bool is_in ();
    void set_pos ();
  };

  /***
   * Adaptive array neigbourhood node iterator
   */

  template<typename aa_type>
  class AANeighbourNodeIterator : public ArrayIterator<typename aa_type::NodeType*, aa_type::Dim>
  {
  public:
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;

  public:
    AANeighbourNodeIterator ();
    AANeighbourNodeIterator (NodeIterator& iit, RectDomType& radius, bool end = false);
    AANeighbourNodeIterator (NodeIterator& iit, RectDomType& radius, int a_idx);
    AANeighbourNodeIterator (const ArrayIterator<typename aa_type::NodeType*, aa_type::Dim>& iit);
  };

  /***
   * Adaptive array neigbourhood node value iterator
   */

  template<typename aa_type>
  class AANeighbourNValueIterator : public ArrayIterator<typename aa_type::NumericType, aa_type::Dim>
  {
  public:
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;

  public:
    AANeighbourNValueIterator ();
    AANeighbourNValueIterator (NodeIterator& iit, RectDomType& center, RectDomType& radius, bool end = false);
    AANeighbourNValueIterator (NodeIterator& iit, RectDomType& center, RectDomType& radius, int a_idx);
    AANeighbourNValueIterator (typename aa_type::NumericArrayType* array, RectDomType& center, 
			       RectDomType& radius, bool end = false);
    AANeighbourNValueIterator (const ArrayIterator<typename aa_type::NumericType, aa_type::Dim>& iit);
  };

  /***
   *
   */
  template <typename aa_type>
  class AANeighbourIterator
  {
  public:
    typedef aa_type AAType;
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;

  protected:
    aa_type* array;
    IndexType center;
    RectDomType radius;
    AANeighbourNodeIterator<aa_type> nit;
    AANeighbourNValueIterator<aa_type> vit;
    bool need_nit;
    
  public:
    AANeighbourIterator () {};
    AANeighbourIterator (NodeIterator& it, IndexType& icenter, RectDomType& iradius, bool local = true);
    AANeighbourIterator (aa_type* array, typename aa_type::cell_index& cell, RectDomType& iradius);
    AANeighbourIterator (NodeIterator& it, IndexType& icenter, RectDomType& iradius, int a_idx, bool local = true);
    AANeighbourIterator (aa_type* array, typename aa_type::cell_index& cell, RectDomType& iradius, int a_idx);
    AANeighbourIterator& operator++ ();
    NumericType& operator* ();

    /// operator++ in one dimension
    AANeighbourIterator& pp ();
    NumericType& center_element ();
    int volume ();
    bool end_p ();
    IndexType local_position ();
    IndexType global_position ();
    IndexType rel_position ();

    AANeighbourNodeIterator<aa_type>& get_nit ();
    AANeighbourNValueIterator<aa_type>& get_vit ();

    NodeIterator node_iterator ();
    NodeIterator center_node_iterator ();

    void next_node ();
    void next_node_pp ();

    void local_reset ();

    aa_type* get_aa ();

    bool clipped_p ();
    bool center_p ();

    void print_dbg ();

  protected:
    bool set_iters (typename aa_type::Node*);
    void update_iters (typename aa_type::Node*);
    void update_iters_pp (typename aa_type::Node*);
    int get_local_pos (typename aa_type::Node*, int center_active_idx, 
		       int& rel_level, int& size, int& center_size, bool swap_coord = false);
    void get_level_size (typename aa_type::Node* element, int active_idx, int& level, int& size);
    IndexType rel_pos (typename aa_type::Node* element);
    void next_node (typename aa_type::Node*);
    void next_node_pp (typename aa_type::Node*);

  };

  template <typename aa_type>
  class AAVirtNeighbourIterator : public AANeighbourIterator<aa_type>
  {
  public:
    typedef aa_type AAType;
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;
    typedef typename aa_type::Node NodeType;

  protected:
    typename aa_type::Node* virtual_node;

  public:
    AAVirtNeighbourIterator (NodeIterator& it, NodeType* ivirual_node, 
			     IndexType& icenter, RectDomType& iradius, bool local = true);
    AAVirtNeighbourIterator (aa_type* array, typename aa_type::cell_index& cell, 
			     NodeType* ivirual_node, RectDomType& iradius);

    AAVirtNeighbourIterator (NodeIterator& it, NodeType* ivirual_node, IndexType& icenter, 
			     RectDomType& iradius, int a_idx, bool local = true);
    AAVirtNeighbourIterator (aa_type* array, typename aa_type::cell_index& cell, 
			     NodeType* ivirual_node, RectDomType& iradius, int a_idx);
    
    AAVirtNeighbourIterator& operator++ ();

    AAVirtNeighbourIterator& pp ();

    void next_node ();
    void next_node_pp ();

    IndexType rel_position ();
    void local_reset ();
  };

  /**
   */
  template <typename input_iterator>
  class AACacheIterator
  {
  public:
    typedef typename input_iterator::AAType::node_iterator NodeIterator;
    typedef typename input_iterator::AAType::NumericType NumericType;
    typedef typename input_iterator::AAType::IndexType IndexType;
    typedef blitz::RectDomain<input_iterator::AAType::Dim> RectDomType;
    typedef typename input_iterator::AAType::Node NodeType;

  protected:
    enum 
      {
	Cnt = 0,
	OffsetOfNode = 1,
	Offsets = 2,
      };

    enum
      {
	ValidCenter,
	InvalidCenter
      };

  protected:
    /**
     * Neigborhood node offset structure:
     *
     * 0. number of offsets for node (size of chunk of offsets)
     * 1. offset of node
     * 2. pointer to offsets of node with refinement 2^0
     * ....
     * AA.refine_level + 2. pointer to offsets of node with refinement 2^(AA.refine_level)
     * 
     */
    
    /**
     * Structure of all refinements of node
     *
     * 0. Offset of node offsets at refinement 2^0
     * 1. Offset of node offsets at refinement 2^1
     * ....
     * AA.refine_level. Offset of node offsets at refinement 2^(AA.refine_level)
     * Offset of node offsets at refinement 2^0 + 0: offset of Neigborhood node offset structure for position 0
     * ....
     * Offset of node offsets at refinement 2^1 + 0: offset of Neigborhood node offset structure for position 0
     * Offset of node offsets at refinement 2^1 + 1: offset of Neigborhood node offset structure for position 1
     * ....
     * Offset of node offsets at refinement 2^(AA.refine_level) + 0: offset of Neigborhood node offset structure for position 0
     * Offset of node offsets at refinement 2^(AA.refine_level) + 1: offset of Neigborhood node offset structure for position 1
     * ....
     * Offset of node offsets at refinement 2^(AA.refine_level) + AA.refine_level - 1: 
     *                              offset of Neigborhood node offset structure for position AA.refine_level - 1
     * ....
     */

    std::vector<std::pair<int, int> > ranges;
    std::vector<int> offsets; 

  public:
    AACacheIterator (typename input_iterator::AAType* array, RectDomType& iradius, 
		     int active_idx, typename input_iterator::AAType::Node* virt_node = NULL);

  protected:
    inline int store_node_offset (AANeighbourNValueIterator<typename input_iterator::AAType>& vit, 
				  int active_idx, std::vector<int>* storage);
    void store_node_offsets (input_iterator& it, std::vector<int>* storage);
    void store_neighborhood_offsets (input_iterator& it, std::vector<int>* storage);
    bool store_offsets (NodeIterator& nit, RectDomType& iradius, int active_idx,
			std::vector<int>* storage, NodeType* virt_node = NULL);
    void make_ranges (typename input_iterator::AAType* aa,  RectDomType& iradius, 
		      int active_idx, NodeType* virt_node = NULL);
  };
 
  template <typename aa_type, typename input_iterator>
  class input_iterator_w
  {
  public:
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;
    typedef typename aa_type::Node NodeType;

  public:
    static input_iterator create (NodeIterator& it, IndexType& icenter, 
				  RectDomType& iradius, int a_idx, NodeType* ivirual_node);
  };

  template <typename aa_type> 
  class input_iterator_w <aa_type, AANeighbourIterator <aa_type> >
  {
  public:
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;
    typedef typename aa_type::Node NodeType;
    
  public:
    static AANeighbourIterator<aa_type> create (NodeIterator& it, IndexType& icenter, 
						RectDomType& iradius, int a_idx, NodeType* ivirual_node)
    { return AANeighbourIterator<aa_type> (it, icenter, iradius, a_idx); }
  };

  template <typename aa_type> 
  class input_iterator_w <aa_type, AAVirtNeighbourIterator<aa_type> >
  {
  public:
    typedef typename aa_type::node_iterator NodeIterator;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::IndexType IndexType;
    typedef blitz::RectDomain<aa_type::Dim> RectDomType;
    typedef typename aa_type::Node NodeType;
    
  public:
    static AAVirtNeighbourIterator<aa_type> create (NodeIterator& it, IndexType& icenter, 
						    RectDomType& iradius, int a_idx, NodeType* ivirual_node)
    { return AAVirtNeighbourIterator<aa_type> (it, ivirual_node, icenter, iradius, a_idx); }
  };


  // *** Implementation of ArrayIterator ***

  template<typename type, int dim>
  ArrayIterator<type, dim>::ArrayIterator
  ()
    :  array (NULL), /* pos, */ sign (true), active_idx (0), clipped (false)
  {
    center = blitz::RectDomain<dim> (blitz::TinyVector<int, dim> (0), blitz::TinyVector<int, dim> (0));
    radius = blitz::RectDomain<dim> (blitz::TinyVector<int, dim> (0), blitz::TinyVector<int, dim> (0));
  }

  template<typename type, int dim>
  ArrayIterator<type, dim>::ArrayIterator
  (blitz::Array<type, dim>* iarray, 
   blitz::RectDomain<dim>& icenter, 
   blitz::RectDomain<dim>& iradius,
   bool end)
    :  array (iarray), center (icenter), radius (iradius), clipped (false)
  {
    if (end)
      { pos = dbc::RectDomIterator<dim> (); sign = true; active_idx = dim; }
    else
      { 
	active_idx = 0; 
	sign = true;
	set_pos ();

	if (!is_in ())
	  ++(*this);
      }
  }

  template<typename type, int dim>
  ArrayIterator<type, dim>::ArrayIterator
  (blitz::Array<type, dim>* iarray, 
   blitz::RectDomain<dim>& icenter, 
   blitz::RectDomain<dim>& iradius,
   int c_idx)
    : array (iarray), center (icenter), radius (iradius), clipped (false)
  { set_active_idx (c_idx); }

  template<typename type, int dim>
  void
  ArrayIterator<type, dim>::set_active_idx (int c_idx)
  {
    active_idx = c_idx; 
    sign = true;
    set_pos ();

    if (!is_in ())
      pp ();
  }

  template<typename type, int dim>
  void
  ArrayIterator<type, dim>::set_pos ()
  {
    typename ArrayIterator<type, dim>::IndexType size (center.ubound () - center.lbound ());
    typename ArrayIterator<type, dim>::IndexType lb (0);
    typename ArrayIterator<type, dim>::IndexType ub (size);

    if (sign)
      {
	lb [active_idx] = (radius.lbound () [active_idx] > 0) ? 0 : radius.lbound () [active_idx];
	ub [active_idx] = (radius.ubound () [active_idx] < 0) ? radius.ubound () [active_idx] : 0;
      }
    else
      {
	lb [active_idx] = (radius.lbound () [active_idx] < 0) 
	  ? size [active_idx] : (radius.lbound () [active_idx] + size [active_idx]);
	ub [active_idx] = (radius.ubound () [active_idx] < 0) 
	  ? size [active_idx] : (radius.ubound () [active_idx] + size [active_idx]);
      }

    blitz::TinyVector<int,dim> tmp (radius.ubound () + (sign ? 0 : size));
    blitz::RectDomain<dim> tmp1 (rect_dom_intersection (radius.lbound (), tmp, lb, ub));

    pos = typename ArrayIterator<type, dim>::RectDomIter (tmp1);
  }

  template<typename type, int dim>
  ArrayIterator<type, dim>::ArrayIterator
  (const ArrayIterator<type, dim>& it)
    : array (it.array), center (it.center), pos (it.pos), sign (it.sign), active_idx (it.active_idx), radius (it.radius)
  {}

  template<typename type, int dim>
  ArrayIterator<type, dim>&
  ArrayIterator<type, dim>::operator++ ()
  {
    while (1)
      {
	++ pos;
	if (pos.end_p () || pos.empty_p ())
	  {
	    if (sign)
	      sign = false;
	    else
	      {
		sign = true;
		active_idx ++;
	      }

	    if (active_idx < dim)
	      set_pos ();
	    else
	      {
		pos = dbc::RectDomIterator<dim> ();
		break;
	      }
	  }

	if (is_in ())
	  break;
      }

    return *this;
  }

  template<typename type, int dim>
  ArrayIterator<type, dim>&
  ArrayIterator<type, dim>::pp ()
  {
    while (1)
      {
	++ pos;
	if (pos.end_p () || pos.empty_p ())
	  {
	    if (sign)
	      {
		sign = false;
		set_pos ();
	      }
	    else
	      {
		pos = dbc::RectDomIterator<dim> ();
		break;
	      }

	    if (is_in ())
	      break;
	  }
      }

    return *this;
  }

  template<typename type, int dim>
  ArrayIterator<type, dim>
  ArrayIterator<type, dim>::operator++ (int)
  {
    ArrayIterator<type, dim> tmp = *this;
    ++(*this);
    return tmp;
  }

  template<typename type, int dim>
  type&
  ArrayIterator<type, dim>::operator* ()
  { return (*array) (blitz::TinyVector<int,dim> (center.lbound () + *pos)); }

  template<typename type, int dim>
  bool
  ArrayIterator<type, dim>::operator== 
  (const ArrayIterator<type, dim>& itt)
  {
    if (array != itt.array)
      return false;

    if (neq_rect_dom<dim> (center, itt.center))
      return false;
    
    if (pos != itt.pos)
      return false;

    if (sign != itt.sign)
      return false;

    if (active_idx != itt.active_idx)
      return false;

    return true;
  }

  template<typename type, int dim>
  bool
  ArrayIterator<type, dim>::operator!= 
  (const ArrayIterator<type, dim>& itt)
  {
    if (array != itt.array)
      return true;

    if (neq_rect_dom<dim> (center, itt.center))
      return true;

    if (pos != itt.pos)
      return true;

    if (sign != itt.sign)
      return true;

    if (active_idx != itt.active_idx)
      return true;

    return false;
  }

  template<typename type, int dim>
  typename ArrayIterator<type, dim>::IndexType&
  ArrayIterator<type, dim>::rel_position ()
  { return *pos; }

  template<typename type, int dim>
  typename ArrayIterator<type, dim>::IndexType
  ArrayIterator<type, dim>::position ()
  { return center.lbound () + *pos; }

  template<typename type, int dim>
  typename ArrayIterator<type, dim>::IndexType&
  ArrayIterator<type, dim>::center_position ()
  { return center.lbound (); }

  template<typename type, int dim>
  type&
  ArrayIterator<type, dim>::center_element ()
  { return (*array) (center.lbound ()); }

  template<typename type, int dim>
  int
  ArrayIterator<type, dim>::active_dim ()
  { return active_idx; }

  template<typename type, int dim>
  bool
  ArrayIterator<type, dim>::valid_p ()
  { return (array != NULL); }

  template<typename type, int dim>
  void
  ArrayIterator<type, dim>::invalidate ()
  { array = NULL; }

  template<typename type, int dim>
  typename ArrayIterator<type, dim>::RectDomType
  ArrayIterator<type, dim>::index_boundary ()
  {
    typedef typename ArrayIterator<type, dim>::IndexType IndexType;

    IndexType lb (center + radius.lbound ());
    IndexType ub (center + radius.ubound ());

    return typename ArrayIterator<type, dim>::RectDomType (lb, ub);
  }

  template<typename type, int dim>
  blitz::Array<type, dim>*
  ArrayIterator<type, dim>::data ()
  { return array; }

  template<typename type, int dim>
  inline bool
  ArrayIterator<type, dim>::end_p ()
  { 
    if (!pos.empty_p ())
      return false; 

    if (sign != true)
      return false;

    if (active_idx != dim)
      return false;

    return true;
  }

  template<typename type, int dim>
  inline bool
  ArrayIterator<type, dim>::clipped_p ()
  { return clipped; }

  template<typename type, int dim>
  bool
  ArrayIterator<type, dim>::is_in ()
  { 
    if (pos.empty_p ())
      return false;

    for (int i = 0; i < dim; ++i)
      {
	int idx = center.lbound () [i] + (*pos) [i];
    
	if ((idx < 0) || (idx >= array->extent (i)))
	  {
	    clipped = true;
	    return false;
	  }
      }

    return true;
  }

  template<typename type, int dim>
  void
  ArrayIterator<type, dim>::print_dbg ()
  {
    std::cout << "center: ";
    print_dbg (center);
    std::cout << "\nposition: ";
    pos.print_dbg ();
    std::cout << "\nradius: ";
    print_dbg (radius);
    std::cout << "sign: " << sign << " active idx: " << active_idx << " clipped: " << clipped << "\n";
  }

  // *** Implementation of AANeighbourNodeIterator ***

  template<typename aa_type>
  AANeighbourNodeIterator<aa_type>::AANeighbourNodeIterator ()
    : ArrayIterator<typename aa_type::NodeType*, aa_type::Dim> () 
  {}

  template<typename aa_type>
  AANeighbourNodeIterator<aa_type>::AANeighbourNodeIterator
  (typename AANeighbourNodeIterator<aa_type>::NodeIterator& iit, 
   typename AANeighbourNodeIterator<aa_type>::RectDomType& radius, bool end)
  {
    typename AANeighbourNodeIterator<aa_type>::RectDomType dom (iit.position (), iit.position () + 1);

    *this = ArrayIterator<typename aa_type::NodeType*, aa_type::Dim> 
      (iit.base (), dom, radius, end);
  }

  template<typename aa_type>
  AANeighbourNodeIterator<aa_type>::AANeighbourNodeIterator
  (typename AANeighbourNodeIterator<aa_type>::NodeIterator& iit, 
   typename AANeighbourNodeIterator<aa_type>::RectDomType& radius, int a_idx)
  {
    typename AANeighbourNodeIterator<aa_type>::RectDomType dom (iit.position (), iit.position () + 1);

    *this = ArrayIterator<typename aa_type::NodeType*, aa_type::Dim> 
      (iit.base (), dom, radius, a_idx);
  }

  template<typename aa_type>
  AANeighbourNodeIterator<aa_type>::AANeighbourNodeIterator
  (const ArrayIterator<typename aa_type::NodeType*, aa_type::Dim>& iit)
    : ArrayIterator<typename aa_type::NodeType*, aa_type::Dim> (iit) 
  {}


  // *** Implementation of AANeighbourNValueIterator ***
  
  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>::AANeighbourNValueIterator ()
    : ArrayIterator<typename aa_type::NumericType, aa_type::Dim> () 
  {}

  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>::AANeighbourNValueIterator
  (typename aa_type::node_iterator& iit, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& center, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& radius,
   bool end)
    : ArrayIterator<typename aa_type::NumericType, aa_type::Dim> (*iit, center, radius, end) 
  {}

  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>::AANeighbourNValueIterator
  (typename aa_type::node_iterator& iit, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& center, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& radius,
   int a_idx)
    : ArrayIterator<typename aa_type::NumericType, aa_type::Dim> (*iit, center, radius, a_idx) 
  {}

  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>::AANeighbourNValueIterator
  (typename aa_type::NumericArrayType* array, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& center, 
   typename AANeighbourNValueIterator<aa_type>::RectDomType& radius,
   bool end)
    : ArrayIterator<typename aa_type::NumericType, aa_type::Dim> (array, center, radius, end) 
  {}

  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>::AANeighbourNValueIterator
  (const ArrayIterator<typename aa_type::NumericType, aa_type::Dim>& iit)
    : ArrayIterator<typename aa_type::NumericType, aa_type::Dim> (iit) 
  {}


  // *** Implementation of AANeighbourIterator ***

  template<typename aa_type>
  AANeighbourIterator<aa_type>::AANeighbourIterator
  (typename aa_type::node_iterator& it, 
   typename AANeighbourIterator<aa_type>::IndexType& icenter, 
   typename AANeighbourIterator<aa_type>::RectDomType& iradius,
   bool local)
    : array (it.get_aa ()), radius (iradius), 
      need_nit (false)
  {
    blitz::RectDomain<aa_type::Dim> dom ((iradius.lbound() >> it.level ()),
					 lshift_ub<int, aa_type::Dim> (iradius.ubound (), it.level ()));

    nit = AANeighbourNodeIterator<aa_type> (it, dom);

    if (local)
      center = icenter;
    else
      center = icenter % (1 << array->refine_level ());

    if (*it)
      {
	dom = blitz::RectDomain<aa_type::Dim> (center, center + 1);
	vit = AANeighbourNValueIterator<aa_type> (it, dom, iradius);
      }
    else
      {
	vit.invalidate ();
	++(*this);
      }
  }

  template<typename aa_type>
  AANeighbourIterator<aa_type>::AANeighbourIterator
  (typename aa_type::node_iterator& it, 
   typename AANeighbourIterator<aa_type>::IndexType& icenter, 
   typename AANeighbourIterator<aa_type>::RectDomType& iradius,
   int a_idx,
   bool local)
    : array (it.get_aa ()), radius (iradius), 
      need_nit (false)
  {
    blitz::RectDomain<aa_type::Dim> dom ((iradius.lbound() >> it.level ()),
					 lshift_ub<int, aa_type::Dim> (iradius.ubound (), it.level ()));

    nit = AANeighbourNodeIterator<aa_type> (it, dom, a_idx);

    if (local)
      center = icenter;
    else
      center = icenter % (1 << array->refine_level ());

    if (*it)
      {
	dom = blitz::RectDomain<aa_type::Dim> (center, center + 1);
	vit = AANeighbourNValueIterator<aa_type> (it, dom, iradius, a_idx);
      }
    else
      {
	vit.invalidate ();
	pp ();
      }
  }

  template<typename aa_type>
  AANeighbourIterator<aa_type>::AANeighbourIterator
  (aa_type* array,
   typename aa_type::cell_index& cell, 
   typename AANeighbourIterator<aa_type>::RectDomType& iradius)
  { 
    typename aa_type::node_iterator it (array->at (cell.node (), false));

    *this = AANeighbourIterator (it, cell.val (), iradius); 
  }

  template<typename aa_type>
  AANeighbourIterator<aa_type>::AANeighbourIterator
  (aa_type* array,
   typename aa_type::cell_index& cell, 
   typename AANeighbourIterator<aa_type>::RectDomType& iradius,
   int a_idx)
  { 
    typename aa_type::node_iterator it (array->at (cell.node (), false));

    *this = AANeighbourIterator (it, cell.val (), iradius, a_idx); 
  }

  template<typename aa_type>
  AANeighbourIterator<aa_type>&
  AANeighbourIterator<aa_type>::operator++ ()
  {
    update_iters (nit.center_element ());
    return *this;
  }

  template<typename aa_type>
  AANeighbourIterator<aa_type>&
  AANeighbourIterator<aa_type>::pp ()
  {
    update_iters_pp (nit.center_element ());
    return *this;
  }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::update_iters (typename aa_type::Node* center_element)
  {
    do {
      if (vit.valid_p ())
	{
	  ++vit;
	  if (! vit.end_p ())
	    return;
	}
      
      if (need_nit)
	++nit;
      else
	need_nit = true;
      
      if (nit.end_p ())
	return;

    } while (!set_iters (center_element));
  }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::update_iters_pp (typename aa_type::Node* center_element)
  {
    do {
      if (vit.valid_p ())
	{
	  if (need_nit)
	    ++vit;
	  else
	    vit.pp ();
	  if (! vit.end_p ())
	    return;
	}
      
      if (need_nit)
	nit.pp ();
      else
	need_nit = true;
      
      if (nit.end_p ())
	return;

    } while (!set_iters (center_element));
  }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::next_node (typename aa_type::Node* center_element)
  {
    do {
      if (need_nit)
	++nit;
      else
	need_nit = true;
      
      if (nit.end_p ())
	return;

    } while (!set_iters (center_element));
  }


  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::next_node ()
  { next_node (nit.center_element ()); }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::next_node_pp (typename aa_type::Node* center_element)
  {
    do {
      if (need_nit)
	nit.pp ();
      else
	need_nit = true;
      
      if (nit.end_p ())
	return;

    } while (!set_iters (center_element));
  }


  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::next_node_pp ()
  { next_node_pp (nit.center_element ()); }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::get_level_size (typename aa_type::Node* element, int active_idx, int& level, int& size)
  {
    if (!element)
      {
	level = 0;
	size = 1;
      }
    else
      {
	level = element->level;
	size = element->data.extent (active_idx);
      }
  }

  template<typename aa_type>
  int
  AANeighbourIterator<aa_type>::get_local_pos 
  (typename aa_type::Node* center_element, int center_active_idx, 
   int& rel_level, int& size, int& center_size, bool swap_coord)
  {
    int nit_level, nit_center_level;
    int active_idx = nit.active_dim ();
    
    get_level_size (*nit, active_idx, nit_level, size);
    get_level_size (center_element, active_idx, nit_center_level, center_size);

    if (swap_coord)
      {
	std::swap (nit_level, nit_center_level);
	std::swap (size, center_size);
      }

    rel_level = nit_level - nit_center_level;
    int local_pos = nit.rel_position () [active_idx] * size * (swap_coord ? -1 : 1);

    if (local_pos < 0)
      {
	if (rel_level < 0)
	  local_pos = -local_pos + (center_active_idx >> (-rel_level));
	else if (rel_level > 0)
	  local_pos = -local_pos + (center_active_idx << rel_level);
	else
	  local_pos = -local_pos + center_active_idx;
      }
    else
      {
	int rel_pos = center_size - center_active_idx;

	if (rel_level < 0)
	  local_pos = -local_pos + size - ((rel_pos >> (-rel_level)) + PMOD1 (rel_pos, -rel_level));
	else if (rel_level > 0)
	  local_pos = -local_pos + size - (rel_pos << rel_level);
	else
	  local_pos = -local_pos + size - rel_pos;
      }

    return local_pos;
  }

  template<typename aa_type>
  bool
  AANeighbourIterator<aa_type>::set_iters (typename aa_type::Node* center_element)
  {
    int rel_level, size, center_size;
    int active_idx = nit.active_dim ();

    int local_pos = get_local_pos (center_element, center [active_idx], rel_level, size, center_size);

    IndexType c;
    IndexType lb (0), ub (1);
    int rlb, rub;
    
    if (rel_level < 0)
      {
	c = center >> (-rel_level);
	rlb = radius.lbound () [active_idx];
	rub = radius.ubound () [active_idx];
      }
    else if (rel_level > 0)
      {
	c = center << rel_level;
	ub = 1 << rel_level;
	rlb = radius.lbound () [active_idx] << rel_level;
	rub = (radius.ubound () [active_idx] << rel_level);
      }
    else
      {
	c = center;
	rlb = radius.lbound () [active_idx];
	rub = radius.ubound () [active_idx];
      }
    
    if (local_pos > 0)
      {
	c [active_idx] = local_pos;
	lb [active_idx] = -local_pos;
	ub [active_idx] = -local_pos + size;
      }
    else
      {
	c [active_idx] = local_pos;
	lb [active_idx] = -local_pos - ((rel_level > 0) ? (1 << rel_level) : 1);
	ub [active_idx] = -local_pos + size - ((rel_level > 0) ? (1 << rel_level) : 1);
      }
    
    if (rel_level < 0)
      {
	int lbr = lb [active_idx] << -rel_level;
	int ubr = ub [active_idx] << -rel_level;

	if (local_pos < 0)
	  {
	    lbr += center_size - center [active_idx] - 1;
	    ubr += center_size - center [active_idx] - 1;
	  }
	else
	  {
	    lbr += center [active_idx];
	    ubr += center [active_idx];
	  }

	if (ubr < rlb)
	  {
	    ub [active_idx] = (rlb >> (-rel_level));
	    lb [active_idx] = (rlb >> (-rel_level));
	  }
	else if (lbr > rub)
	  {
	    lb [active_idx] = (rub >> (-rel_level));
	    ub [active_idx] = (rub >> (-rel_level));
	  }
      }
    else
      {
	if (ub [active_idx] < rlb)
	  {
	    ub [active_idx] = rlb;
	    lb [active_idx] = rlb;
	  }
	else if (lb [active_idx] > rub)
	  {
	    lb [active_idx] = rub;
	    ub [active_idx] = rub;
	  }
	else
	  {
	    if (ub [active_idx] > rub)
	      ub [active_idx] = rub;
	    if (lb [active_idx] < rlb)
	      lb [active_idx] = rlb;
	  }
      }

    RectDomType radius (lb, ub);
    RectDomType cr (c, c + ((rel_level > 0) ? (1 << rel_level) : 1));

    if (rect_dom_empty_p (radius))
      {
	vit.invalidate ();
	return false;
      }

    if (!*nit)
      vit.invalidate ();
    else
      {
	vit = AANeighbourNValueIterator<aa_type> (&(*nit)->data, cr, radius);
	if (vit.end_p ())
	  return false;
      }

    return true;
  }

  template<typename aa_type>
  typename aa_type::NumericType&
  AANeighbourIterator<aa_type>::operator* ()
  {
    if (!vit.valid_p ())
      {
	typename aa_type::IndexType tmp (nit.position ());
	return array->data (tmp);
      }
    return *vit;
  }

  template<typename aa_type>
  typename aa_type::NumericType&
  AANeighbourIterator<aa_type>::center_element ()
  {
    if (!vit.valid_p ())
      {
	typename aa_type::IndexType tmp (nit.center_position ());
	return array->data (tmp);
      }
    return vit.center_element ();
  }

  template<typename aa_type>
  int
  AANeighbourIterator<aa_type>::volume ()
  {
    if (!need_nit)
      {
	if (!nit.center_element ())
	  return power (1 << array->refine_level (), aa_type::Dim);
	return power (1 << (array->refine_level () - (nit.center_element ())->level), aa_type::Dim);
      }

    if (!*nit)
      return power (1 << array->refine_level (), aa_type::Dim);
    return power (1 << (array->refine_level () - (*nit)->level), aa_type::Dim);
  }

  template<typename aa_type>
  bool
  AANeighbourIterator<aa_type>::end_p ()
  { return nit.end_p () && (!vit.valid_p () || (vit.valid_p () && vit.end_p ())); }

  template<typename aa_type>
  typename aa_type::IndexType
  AANeighbourIterator<aa_type>::local_position ()
  { return vit.valid_p () ? vit.position () : typename aa_type::IndexType (0); }

  template<typename aa_type>
  typename aa_type::IndexType
  AANeighbourIterator<aa_type>::global_position ()
  { return need_nit ? nit.position () : nit.center_position (); }

  template<typename aa_type>
  typename aa_type::IndexType
  AANeighbourIterator<aa_type>::rel_position ()
  { return rel_pos (nit.center_element ()); }

  template<typename aa_type>
  typename aa_type::IndexType
  AANeighbourIterator<aa_type>::rel_pos (typename aa_type::Node* center_element)
  { 
    if (vit.valid_p () && !need_nit)
      return vit.rel_position ();

    int rel_level, size, center_size, local_pos;
    int active_idx = nit.active_dim ();
    typename aa_type::IndexType rel_pos (0);

    if (vit.valid_p ())
      local_pos = get_local_pos (center_element, vit.position () [active_idx],
				rel_level, center_size, size, true);
    else
      local_pos = get_local_pos (center_element, 0,
				 rel_level, center_size, size, true);
    
    if (local_pos > 0)
      rel_pos [active_idx] = local_pos + (center_size - center [active_idx]);
    else
      rel_pos [active_idx] = local_pos + center [active_idx] + size - 1;

    return rel_pos;
  }

  template<typename aa_type>
  AANeighbourNodeIterator<aa_type>&
  AANeighbourIterator<aa_type>::get_nit ()
  { return nit; }

  template<typename aa_type>
  AANeighbourNValueIterator<aa_type>&
  AANeighbourIterator<aa_type>::get_vit ()
  { return vit; }

  template<typename aa_type>
  typename aa_type::node_iterator
  AANeighbourIterator<aa_type>::node_iterator ()
  { 
    typename aa_type::IndexType tmp (need_nit ? nit.position () : nit.center_position ());
    return array->at (tmp, false); 
  }

  template<typename aa_type>
  typename aa_type::node_iterator
  AANeighbourIterator<aa_type>::center_node_iterator ()
  { return array->at (nit.center_position (), false); }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::local_reset ()
  { set_iters (nit.center_element ()); }

  template<typename aa_type>
  aa_type*
  AANeighbourIterator<aa_type>::get_aa ()
  { return array; }

  template<typename aa_type>
  bool
  AANeighbourIterator<aa_type>::clipped_p ()
  { return nit.clipped_p (); }

  template<typename aa_type>
  bool
  AANeighbourIterator<aa_type>::center_p ()
  { return nit.center_element () && vit.valid_p () ; }

  template<typename aa_type>
  void
  AANeighbourIterator<aa_type>::print_dbg ()
  {
    std::cout << "center: ";
    print_dbg (center);
    std::cout << "\nradius: ";
    print_dbg (radius);
    std::cout << "\nnit:";
    nit.print_dbg ();
    std::cout << "\nvit:";
    vit.print_dbg ();
    std::cout << "\nneed_nit: " << need_nit << "\n";
  }

  // *** Implementation of AAVirtNeighbourIterator ***

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>::AAVirtNeighbourIterator
  (typename aa_type::node_iterator& it, 
   typename aa_type::Node* ivirtual_node,
   typename AAVirtNeighbourIterator<aa_type>::IndexType& icenter, 
   typename AAVirtNeighbourIterator<aa_type>::RectDomType& iradius,
   bool local)
    : virtual_node (ivirtual_node)
  {
    BOOST_ASSERT (ivirtual_node != NULL);
    BOOST_ASSERT (ivirtual_node->level >= it.level ());

    this->array = it.get_aa (), this->radius = iradius,  this->need_nit = false;

    blitz::RectDomain<aa_type::Dim> dom ((iradius.lbound() >> virtual_node->level),
					 lshift_ub<int, aa_type::Dim> (iradius.ubound (), virtual_node->level));

    this->nit = AANeighbourNodeIterator<aa_type> (it, dom);

    if (local)
      this->center = icenter;
    else
      this->center = icenter % (1 << this->array->refine_level ());

    if (*it)
      {
	typename AAVirtNeighbourIterator<aa_type>::IndexType c (this->center >> (virtual_node->level-it.level ()));
	dom = blitz::RectDomain<aa_type::Dim> (c, c + 1);
	blitz::RectDomain<aa_type::Dim> rad 
	  (iradius.lbound () >> (virtual_node->level - it.level ()),
	   lshift_ub<int, aa_type::Dim> (iradius.ubound (), virtual_node->level - it.level ()));

	this->vit = AANeighbourNValueIterator<aa_type> (it, dom, rad);
      }
    else
      {
	this->vit.invalidate ();
	++(*this);
      }
  }

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>::AAVirtNeighbourIterator
  (typename aa_type::node_iterator& it, 
   typename aa_type::Node* ivirtual_node,
   typename AAVirtNeighbourIterator<aa_type>::IndexType& icenter, 
   typename AAVirtNeighbourIterator<aa_type>::RectDomType& iradius,
   int a_idx,
   bool local)
    : virtual_node (ivirtual_node)
  {
    BOOST_ASSERT (ivirtual_node != NULL);
    BOOST_ASSERT (ivirtual_node->level >= it.level ());

    this->array = it.get_aa (), this->radius = iradius,  this->need_nit = false;

    blitz::RectDomain<aa_type::Dim> dom ((iradius.lbound() >> virtual_node->level),
					 lshift_ub<int, aa_type::Dim> (iradius.ubound (), virtual_node->level));

    this->nit = AANeighbourNodeIterator<aa_type> (it, dom, a_idx);

    if (local)
      this->center = icenter;
    else
      this->center = icenter % (1 << this->array->refine_level ());

    if (*it)
      {
	typename AAVirtNeighbourIterator<aa_type>::IndexType c (this->center >> (virtual_node->level-it.level ()));
	dom = blitz::RectDomain<aa_type::Dim> (c, c + 1);
	blitz::RectDomain<aa_type::Dim> rad 
	  (iradius.lbound () >> (virtual_node->level - it.level ()),
	   lshift_ub<int, aa_type::Dim> (iradius.ubound (), virtual_node->level - it.level ()));

	this->vit = AANeighbourNValueIterator<aa_type> (it, dom, rad, a_idx);
      }
    else
      {
	this->vit.invalidate ();
	pp ();
      }
  }

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>::AAVirtNeighbourIterator
  (aa_type* array,
   typename aa_type::cell_index& cell, 
   typename aa_type::Node* ivirual_node,
   typename AAVirtNeighbourIterator<aa_type>::RectDomType& iradius)
  {
    typename aa_type::node_iterator it (array->at (cell.node (), false));
    
    *this = AAVirtNeighbourIterator (it, ivirual_node, cell.val (), iradius);
  }

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>::AAVirtNeighbourIterator
  (aa_type* array,
   typename aa_type::cell_index& cell, 
   typename aa_type::Node* ivirual_node,
   typename AAVirtNeighbourIterator<aa_type>::RectDomType& iradius,
   int a_idx)
  {
    typename aa_type::node_iterator it (array->at (cell.node (), false));
    
    *this = AAVirtNeighbourIterator (it, ivirual_node, cell.val (), iradius, a_idx);
  }

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>&
  AAVirtNeighbourIterator<aa_type>::operator++ ()
  {
    update_iters (virtual_node);
    return *this;
  }

  template<typename aa_type>
  AAVirtNeighbourIterator<aa_type>&
  AAVirtNeighbourIterator<aa_type>::pp ()
  {
    update_iters_pp (virtual_node);
    return *this;
  }

  template<typename aa_type>
  typename aa_type::IndexType
  AAVirtNeighbourIterator<aa_type>::rel_position ()
  { return rel_pos (virtual_node); }

  template<typename aa_type>
  void
  AAVirtNeighbourIterator<aa_type>::next_node ()
  { next_node (virtual_node); }

  template<typename aa_type>
  void
  AAVirtNeighbourIterator<aa_type>::next_node_pp ()
  { next_node_pp (virtual_node); }

  template<typename aa_type>
  void
  AAVirtNeighbourIterator<aa_type>::local_reset ()
  { set_iters (virtual_node); }

  
  // *** Implementation of AACacheIterator ***

  template<typename input_iterator>
  AACacheIterator<input_iterator>::AACacheIterator
  (typename input_iterator::AAType* array, RectDomType& iradius, int active_idx,
   typename input_iterator::AAType::Node* virt_node)
  { make_ranges (array, iradius, active_idx, virt_node); }

  template<typename input_iterator>
  void
  AACacheIterator<input_iterator>::make_ranges 
  (typename input_iterator::AAType* aa, 
   RectDomType& iradius, 
   int active_idx,
   typename input_iterator::AAType::Node* virt_node)
  {
    typename input_iterator::AAType::node_iterator nit1 (aa);
    
    for (; nit1.end_p (); nit1.pp (active_idx))
      {
	std::pair<int, int> i (nit1.position (active_idx), offsets.size ());
	ranges.push_back (i);
	if (!store_offsets (nit1, iradius, active_idx, &offsets, virt_node))
	  break;
      }

    if (!nit1.end_p ())
      {
	typename input_iterator::AAType::node_iterator nit2 (aa, true);
	std::vector<std::pair<int, int> > tmp;
	
	for (; nit2.end_p (); nit2.mm (active_idx))
	  {
	    std::pair<int, int> i (nit2.position (active_idx), offsets.size ());
	    tmp.push_back (i);
	    if (!store_offsets (nit2, iradius, active_idx, &offsets, virt_node))
	      break;
	  }

	while (tmp.empty ())
	  {
	    ranges.push_back (tmp.back ());
	    tmp.pop_back ();
	  }
      }
  }

  template<typename input_iterator>
  bool
  AACacheIterator<input_iterator>::store_offsets 
  (typename input_iterator::AAType::node_iterator& nit,
   RectDomType& iradius, 
   int active_idx,
   std::vector<int>* storage,
   typename input_iterator::AAType::Node* virt_node)
  {
    typedef typename input_iterator::AAType aa_type;
    typedef blitz::TinyVector<int, aa_type::Dim> IndexType;
    typedef typename aa_type::NumericType NumericType;

    bool clipped = false;
    typename aa_type::Node* bckp = nit.node ();
    typename aa_type::Node tmp_node (1);
    int scale = 1;
    int begining = storage->size ();
    storage->resize (nit.get_aa ()->refine_level ());

    nit.node () = NULL;

    for (int i = 0; i < nit.get_aa ()->refine_level (); ++i, scale <<= 1)
      {
	IndexType pos (0);
	int local_begining = storage->size ();

	(*storage) [begining + i] = local_begining;
	storage->resize (scale);
	
	for (int j = 0; j < scale; ++j)
	  {
	    input_iterator it;

	    (*storage) [local_begining + j] = storage->size ();

	    pos [active_idx] = j;
	    it = input_iterator_w<typename input_iterator::AAType, input_iterator>::create (nit, pos, iradius, active_idx, virt_node);

	    store_neighborhood_offsets (it, storage);

	    if (it.clipped_p ())
	      clipped = true;
	  }

	tmp_node.data.resize (IndexType (scale << 1));
	nit.node () = &tmp_node;
      }

    nit.node () = bckp;

    return clipped;
  }

  template<typename input_iterator>
  inline int
  AACacheIterator<input_iterator>::store_node_offset (AANeighbourNValueIterator<typename input_iterator::AAType>& vit, 
						      int active_idx, std::vector<int>* storage)
  {
    int cnt = 0;

    for (vit.set_active_idx (active_idx);
	 vit.end_p (); vit.pp (), ++cnt)
      {
	int idx = a_offset<typename input_iterator::AAType::NumericType, input_iterator::AAType::Dim> 
	  (*vit.data (), vit.rel_position ());
	storage->push_back (idx);
      }

    return cnt;
  }

  template<typename input_iterator>
  void
  AACacheIterator<input_iterator>::store_node_offsets (input_iterator& it, std::vector<int>* storage)
  {
    typedef typename input_iterator::AAType aa_type;
    typedef blitz::TinyVector<int, aa_type::Dim> IndexType;
    typedef typename aa_type::NumericType NumericType;
    typedef typename aa_type::NodeType NodeType;

    int begining, cnt;
    AANeighbourNValueIterator<aa_type> vit;

    cnt = 0;
    begining = storage->size ();
    storage->resize (begining + OffsetOfNode + it.get_aa ()->refine_level ());
    
    (*storage) [begining + OffsetOfNode] = 
      a_offset<NodeType*, aa_type::Dim> (*(it.get_nit ().data ()), 
					 it.get_nit ().rel_position ());

    typename aa_type::Node* bckp = *(it.get_nit ());
    typename aa_type::Node tmp_node (1);

    *(it.get_nit ()) = NULL;

    for (int i = 0; i < it.get_aa ()->refine_level (); ++i)
      {
	it.local_reset ();

	(*storage) [begining + Offsets + i] = storage->size ();
	cnt += store_node_offset (it.get_vit (), it.get_nit ().active_dim (), storage); 

	tmp_node.data.resize (IndexType (2 << i));
	*(it.get_nit ()) = &tmp_node;
      }

    *(it.get_nit ()) = bckp;

    (*storage) [begining + Cnt] = cnt;
  }

  template<typename input_iterator>
  void
  AACacheIterator<input_iterator>::store_neighborhood_offsets (input_iterator& it, std::vector<int>* storage)
  {
    if (it.center_p ())
      {
	storage->push_back (ValidCenter);
	store_node_offset (it.get_vit (), it.get_nit ().active_dim (), storage);
	it.next_node_pp ();
      }
    else
      storage->push_back (InvalidCenter);
    
    for (; it.end_p (); it.next_node_pp ())
      store_node_offsets (it, storage);
  }

}; // namespace dbc

#undef PMOD1

#endif
